The boundary integral method for magnetic billiards 



Klaus HornbergerjJ and Uzy Smilanskyf 

f Department of Physics of Complex Systems, The Weizmann Institute of 
Science, 76100 Rehovot, Israel 

X Max-Planck-Institut fur Physik komplexer Systeme, Nothnitzer Strafie 38, 
01187 Dresden, Germany 



Abstract. We introduce a boundary integral method for two-dimensional 
quantum billiards subjected to a constant magnetic field. It allows to calculate 
spectra and wave functions, in particular at strong fields and scmiclassical values 
of the magnetic length. The method is presented for interior and exterior problems 
with general boundary conditions. We explain why the magnetic analogues of the 
field-free single and double layer equations exhibit an infinity of spurious solutions 
and how these can be eliminated at the expense of dealing with (hyper-)singular 
operators. The high efficiency of the method is demonstrated by numerical 
calculations in the extreme semiclassical regime. 
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1. Introduction 

Magnetic billiards are systems of a confined, charged particle in a constant magnetic 
field. In mesoscopic physics they serve as models to explain shape-dependent features 
of nano-scale devices |2|, like quantum dots. In quantum chaos they are studied as 
natural extensions of planar billiards [[| ||, |[ |(| . These systems are particularly suited 
for the study of semiclassical effects (both, theoretically f?], |[ |9| and in experiments 
[Tl| , [l^| ) since the field strength which essentially determines the scale of quantum 
effects is a free parameter. 

The presence of a Lorentz force severely affects the classical, two-dimensional 
billiard dynamics. The criteria for hyperbolicity are altered ||. For strong enough 
fields closed cyclotron orbits occur, while other trajectories perform a skipping motion 
along the billiard boundary. Most significantly, the exterior dynamics where the 
billiard boundary acts as an obstacle from outside is not a scattering problem like 
in the field free case but exhibits bounded skipping motion around the billiard. 

The magnetic quantum spectra and wave functions reflect these classical 
properties. For strong fields a separation takes place in the spectrum. Close to 
the energies of the Landau levels one finds bulk states which correspond to a free 
cyclotron motion of the particle. In addition, edge states appear which are localized 
at the boundary, corresponding to a skipping motion along it. Unlike the field free 
case, the spectrum is purely discrete also in the exterior, with accumulation points at 
the energies of the Landau levels. 
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From a technical point of view, calculations of spectra and wave functions are 
considerably more difficult with a magnetic field present. So far, they were mostly 
realized by diagonalizing the Hamiltonian jl3[ [hi], |T^, [T(|. This requires the choice 
and truncation of a basis which is problematic in the general case when no natural 
basis exists. It explains why calculations of exterior wave functions were not even 
attempted. 

The spectra of field free billiards are usually calculated by transforming the 
eigenvalue problem into an integral equation of lower dimension. The corresponding 
integral operator is defined in terms of the free Green function and depends only 
on the boundary. This boundary integral method is known to be more efficient than 
diagonalization by an order of magnitude and avoids the arbitrariness of choosing and 
truncating a basis. 

It seems natural to extend these ideas to the magnetic problem. A step in this 
direction was taken by Tiago et al ]l7j who essentially propose a null-field method Q 
which involves the irregular Green function in the angular momentum decomposition. 
A drawback of their approach is that the latter function must be known for large 
angular momenta which is practically inaccessible numerically. 

In the present paper we propose a boundary integral method for two-dimensional 
magnetic billiards. It involves the regular Green function in the position space 
representation. We derive the method for both, the interior and the exterior problem 
and for general boundary conditions which include the Dirichlet and Neumann choice 
as special cases. The method allows for the first time to calculate spectra and wave 
functions of magnetic billiards for arbitrary fields and semiclassical values of the 
magnetic length. Thousands of consecutive energy levels may be calculated to high 
precision with moderate numerical effort. 

1.1. Outline 

For field-free billiards two independent boundary integral equations are known. In 
section 2 we derive their magnetic analogues in a gauge-invariant formulation. It is 
shown that unlike in the field-free case each of these equations yield only a necessary 
but not a sufficient condition for the definition of the spectra. In other words, each 
equation admits spurious solutions. We will identify the physical origin of the latter 
and propose a way to avoid them at the expense of dealing with singular (and possibly 
even hypersingular [fl9||) operators. 

The explicit form of the integral operators is presented in section 3 where we 
discuss the nature of the singularities, too. In section 4 it is shown how the integral 
equations may be solved treating the singular parts of the operators analytically. 
This leaves the remaining problem in a form suitable for numerical treatment. Its 
implementation is sketched in section || together with a discussion of the numerical 
convergence and the attainable accuracy. 

Finally, we demonstrate the power of the proposed method in section [| where we 
study spectral statistics using several thousand levels and present interior and exterior 
wave functions in the quasi-classical regime. 

1.2. Preliminaries 

We are interested in the spectrum of a charged particle constrained to a compact 
domain T> GM 2 which is subject to a constant, perpendicular magnetic field of strength 
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B. Alternatively, one may consider the complementary situation by constraining the 
particle to the exterior R 2 \ T>. Unlike the non-magnetic case the exterior spectrum is 
discrete, with accumulation points at the Landau levels. The stationary Schrodinger 
equation reads 

-^-iftv r - 9 A(r)) 2 ^(r)=i^(r). (1) 
2m 

where to, q, and E are mass, charge, and energy of the particle, respectively. The 
vector potential may be written in terms of the symmetric gauge A sym , 

A(r) = A sym (r) + v*(r) := | r + v X (r), (2) 

where X accounts for the gauge freedom which we shall limit to Coulomb type 
v 2 X(r)=0. 

We assume the domain boundary T — dT> to be smooth and choose its normals 
n(r) to point outwards (i.e. into R 2 \ T>.) Keeping their orientation fixed will allow to 
distinguish the interior and the exterior problem. 

The number of parameters (E, B, h, q, to) in (||) can be reduced. Scaling time by 
the Larmor frequency u> = qB/(2m), one remains with two length scales, 

p 2 = -^ and 6 2 = A (3 ) 

as the only parameters describing the system, p is the classical cyclotron radius 
whereas the magnetic length b has a pure quantum meaning. It gives the mean radius 
of a bulk ground state. The scaled energy may be expressed in terms of the spacing 
between Landau levels 

The expression for the (unsealed) wave number k — \j2mE jh — 2p/b 2 shows that there 
are two distinct short-wave limits: The high-energy limit p — > oo and the semiclassical 
limit 6^0. The former corresponds to increasing the energy at fixed magnetic 
field while in the semiclassical limit one increases both energy and field, keeping p 
fixed. It is important to distinguish between these limits. The high energy direction 
is the simpler one because the dynamical effect of the magnetic field tends to vanish. 
However, for semiclassical studies the later direction is the only proper choice because 
it leaves the classical dynamics unaffected. 

For most of the numerical demonstrations in section ^ the magnetic length b is 
chosen as the spectral variable in order to present the boundary integral method in 
the nontrivial limit. Therefore, to show clearly the dependence of the equations on 
b we do not introduce dimensionless variables for the scaled positions r/b. However, 
we facilitate the replacement by dimensionless variables by stating all expressions 
(including the scaled gradient v r /i, := 6v r ) in terms of that quotient. 



2. Derivation of the boundary integral equations 

In this section two magnetic boundary integral equations are derived. We show 
why they have spurious solutions and how to avoid this by constructing a combined 
boundary integral equation. 



V.' = ± T (d n/b V-iA n V>), rel\ (6) 
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2.1. Single and double layer eguations 

The quantum wave function ip £ /^(R 2 ) is denned by its differential equation 

±hiv r/6 -A(r)] 2 -2^(r) = 0, (5) 

and a specification of boundary conditions. We shall employ general gauge invariant 
boundary conditions, 

A 
b[ 

with A = 2/(Sfe)A(r), A n = n(r)A, d n / b := n(r)v,./ b := 6n(r)v r . Here, A (which 
may be a function of the position on the boundary) interpolates between Dirichlet 
boundary conditions (A = 0) and the Neumann case (A -1 = 0). Equation (||) is a 
generalization of the mixed boundary conditions known for the Helmholtz problem 
pp[ pl| p2[ . It implies that the normal component of the current density operator 
j n = lm(tfj*d n /i ) tp) — A„|-0| 2 vanishes for any A. The lower sign in (|6|) corresponds to 
the exterior problem. 

We mention in passing that another type of boundary conditions for magnetic 
billiards was proposed recently |2j| . 

The Green function satisfies the inhomogeneous equation 

i[-iv r/6 - A(r)] 2 - 21/) G(r; r ) = -\ S(^) . (7) 

Its properties are described in the appendix. Note, that it does not depend on the 
difference vector (r — ro) alone but has a gauge dependent phase, 

G(r;r )=e- i (f 1 -*!*o)) G°(^). (8) 

We take G to be the free regular Green function by demanding 

lim G° v (z) = 0, (9) 

z — >oo 

which specifics G uniquely as the Fourier transform of the free time evolution operator. 
As one expects, the regular Green function decays exponentially once the points are 
separated by a distance, |(r — r )| > 2p, which cannot be traversed classically. An 
independent solution to (Q) exists which grows exponentially beyond this classically 
allowed region. It may be called irregular free Green function and was used in the 
null-field method approach for reasons to be explained below. In the following 
only the regular Green function will be used. 

We start by considering the interior problem. The treatment of the exterior case 
is sketched afterwards. Equations (||) and (0) can be combined to yield a form suitable 
for the Green and Gaufi integral theorems, 

^v 2 /b G-Gv r 2 /b V + 2iv r/b (A^G) ^iPS^) , (10) 

where the bar indicates complex conjugation. Choosing ro £ M 2 \ V the integral of 
( [To|) over the domain T> may be transformed to a boundary integral, 

[i,d n/b G-Gd n/b i> + 2iKi>G}^- = { ^ $21%^ (ID 

corresponding to the interior problem. Next, the vector potential part of the integrand 
is split, 

[V,(a n/i ,G + iA„G)-G(a n/ ^-iA„^)]^ = { ^ ( Q ro) 5^oGR 2 \Z>. ( 12 ) 
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which will allow for a gauge invariant formulation of the boundary integral equation. 
Taking r G T, := r ± en with e > 0, we add the two equations above to obtain 

[i, (cf n/b G + i A„ G £ ) - G £ (d n/b ^ - i A„ 0))] — = |V(ro )• (13) 

Here we have introduced the abbreviations G^ = |G(r;rg) + |G(r;rg), d n , b G = 

\d n / b G(r; Tq) + |9 n /&G(r; rjT ). Equation ( |l3| ) is true for all (sufficiently small) e > 
which means that the limit e — > exists. Moreover, it can be shown that one is allowed 
to exchange the integration with the limit (G — > G, 9, 6 G — * 9„/jG.) Observing the 

boundary condition (|^) and renaming the unknown function, u — d n / b ip — iA n ip, 
u a := u(r ) we get 

[G- ^(9 n/fe G + iA n G)] u^- = ^(-|«o), (14) 

which is an integral equation defined on the boundary T. 

In order to obtain the corresponding equation for the exterior problem consider a 
large disk K v D T> of radius p and integrate (|l(]) over KC\T>. Once r is in the vicinity 
of r, the contribution of dK, to the boundary integral vanishes as p — > oo due to the 
exponential decay of G (since ip € £2)- Similar to ( [l3|) one obtains an equation 

[^{W n/b G + i A„ G £ ) - G £ (9„ /b V - i A„ */,))] ^ = i^(r+) (15) 

which allows for the limit e — > to be taken before performing the integration. The 
resulting boundary integral equation differs from ( |l4| ) only by a sign. In the following, 
we treat both cases simultaneously with the convention that the upper sign stands for 
the interior problem and the lower sign for the exterior case. 

J [GT ^(d„ /b G + iA„G)] u^- = ± (-|«o). (16) 

For historical reasons pif , we will refer to these equations as the single layer equations 
for the interior and the exterior domain. 

A second kind of boundary integral equations can be derived by applying the 
differential operator (d no / b — iA„ ) := n(r )(v ro /f ) — iA(r ) on equations ( |l3| ) and 

TBI), 

V ) (9„ o/6 -iA„ n )(a„ /b G + iA„G )— - J (d no/b G-iA no G )(d n/b ip - ik n iji) — 

= ±i(a„ o/b -iA no )V(4). (17) 

This equation is true for all e > which means that the limit e — > exists. As for 
the first integral, we are again allowed to permute the limit and the integration which 
yields a proper integral. Consequently, the limit of the second integral is finite, too. 
However, in the second integral we are not allowed to exchange the integration with 
taking the limit because the limiting integrand has a l/(r — ro) 2 -singularity which is 
not integrable. 

Integral operators of this kind are named hyper singular []l9| . Similar to a Cauchy 
principal value integral, they are defined by taking a special limit. However, in the 
present case the singularity is stronger by one order. In the next section, we define 
which limit is to be taken. It is denoted by and should be read "finite part of the 
integral" . With this concept and @ we obtain the double layer equations, 

J^(d no/b G-iA no G)u— T ^^(9„ o/6 -iA, i0 )(a il/b G + iA„G)w— = t\u (18) 
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which are again integral equations defined on the boundary F. 

It is useful to introduce a set of integral operators, (whose labels D and N indicate 
correspondence to pure Dirichlet or Neumann conditions) 



drG 



u. 



dr,„ ^ 



Qsi M - / - T (d n/b G + iA n G)u, 



(19) 



dr 



(d no/b G-iA na G)u, Q N 



This way, the requirement of the existence of nontrivial solutions of equations ( |l6|) 
and ([l8]) is equivalent to demanding that the corresponding Fredholm determinants 
vanish. 



dr 



{pn /b ~ iA„ )((9 n / & G + i A„ G)u. 



det 



det 



QslTAQS 



q2 t aq£ ± iid 



(single layer,) 



(double layer. 



(20) 



(21) 



These are secular equations although the explicit dependence on the spectral variable 
is not shown in our abbreviated notation. If one chooses p as the spectral variable, 
only the energy parameter v = p 2 /b 2 of the Green function is varied. Taking b as 
spectral variable will in addition change the intrinsic length scale. 

As mentioned already, each of the determinants ( ^0| ) and (|2l]) may have zeros 
which do not correspond to solutions of the original eigenvalue problem given by ((bJ) 
and (||). For finite e the equations jl^), ([l5|), and ( |l7| ) are still equivalent to the latter. 
They acquire additional spurious solutions only as they are transformed to boundary 
integral equations by the limit e — > 0. 



2.2. Spurious solutions and the combined operator 

The physical origin of the redundant zeros is apparent in our gauge invariant 
formulation. They are proper solutions for the domain complementary to the one 
considered. This is obvious for the single layer equation with Dirichlet boundary 
conditions (A = 0) where the spectral determinant does not depend on the orientation 
of the normals. The same is true for the double layer equation with Neumann boundary 
conditions (A -1 = 0). 

In general, the character of the spurious solutions may be summarized as follows: 
Independently of the boundary conditions, the single layer equation includes the 
Dirichlet solutions of that domain which is complementary to the one considered. 
Likewise, the double layer equation is polluted by the Neumann solutions of the 
complementary domain, irrespective of the boundary conditions employed. 

The statement is easily proved by observing that the single-layer-Neumann 
operator and the double-layer-Dirichlet operator are adjoint to each other, 
Q^} = (Qdi)^ while the operators and are self-adjoint. This is shown explicitly 
in the next section. Now assume that it is a complementary Dirichlet solution. In 
Dirac notation, 

Qg|«)=0 A Q»\u) T 1 1 \u)=i) (22) 

=>■ HQS = ° a («|QS=FiH = o. 

Applying the dual of u to the single layer operator yields 

( U |QSta{hQN T Ih}=0, (23) 
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which implies that the Fredholm determinant of the single layer operator vanishes. 
Similarly, if u is a complementary Neumann solution, 

±Q» + i|u> = A Q»=0 (24) 

=> ±hq£ + ±h = o a («|qS = o 

then its dual satisfies the double layer equation, again for any A, 

±{±HQ£+±H} T AHQ3i = 0. (25) 

Since the spurious solutions are never of the same type it is possible to dispose of them 
by requiring that both the single and the double layer equations should be satisfied 
by the same solution u. Therefore, one obtains a necessary and sufficient condition 
for the definition of the spectrum by considering a combined operator 

Qf := (OS T AQS ± hd) + ia fcfi T A<$ + ^id) . (26) 



2 J \ 81 1 31 2, 

with an arbitrary constant a. It has a zero eigenvalue only if both, single and double 
layer operators do so. The choice ( ^6|) works very well in practice as will be shown 
below. 

It seems natural to require that both the single and the double layer equation must 
be satisfied to determine a proper eigenvalue. The original equation ( pd| ) consists of 
two independent conditions (ro € V and ro £ I 2 \ T>.) Only for special situations, 
such as the field free problem, the two conditions are equivalent so that each of them 
singly provides the correct spectrum. For a discussion of the field free case see e.g. 

It is interesting to note that (for the interior problem) spurious solutions will 
not appear if one uses the irregular Green function. The reason is that the gauge- 
independent part of this function is complex which destroys the mutual adjointness 
of the operators. This is why the irregular Green function had to be chosen for the 
null-field method employed in JlTj]. For the boundary integral method the option to 
use this exponentially divergent solution of (Q) is excluded since the corresponding 
operator would get arbitrarily ill-conditioned once the size of the boundary exceeds 
the cyclotron diameter. 

Our last remark is concerned with the important case of Dirichlet boundary 
conditions. Here, one could as well derive a pair of boundary integral equations that 
are not gauge-invariant. (Just set ip = in (|ll| ) and consider u = d n /bif).) Of course, 
these equations would yield the proper gauge-invariant eigen-energies of the problem. 
However, the energies of the additional spurious solutions would depend on the chosen 
gauge and a characterization of the latter in terms of solutions of a complementary 
problem would not be possible. 

3. The boundary integral operators 

In this section we give explicit expressions for the boundary integrals. This allows to 
define the "finite part integral" appearing in the double layer equation. 

3.1. Explicit expression for the integral kernels 

The form of the Green function (ph leads to the following expressions for the integral 
kernels Q(r;r ) of the operators (|l9|), (Q[u])(r ) = J r drQ(r; r Q )u(r), with n = n(r), 



Boundary integral method for magnetic billiards 



n = n(r ), Ar := (r — r ), and z :— Ar 2 jb 2 



Qf l (r;r ) = 
QS(r;r ) = 

QS(r;r ) = 

Qdi(r;r ) = 



c i(^-XM+X(ro)) G (z) 



jf^-XW+jffro)) 



i^G«(z) + 2^zi : G0 W } > 



6 2 

. Ar x no 



dz 



'rxro 



-X( r )+X( r o), 



6 2 Ar 2 dz 

(Ar x n )(Ar x n) . n x n 



b 4 



b 2 



+(2i 



. n x no 



b 2 



Ar z dz 



i (Arn)(Arn ) ^ 



G°(z) 
d 2 



(27) 
(28) 

(29) 

(30) 



rG 



Ar* 



Note, that the gauge freedom X cancels in the prefactors and only appears in the phase. 
It can be absorbed by the function u(r) — ► exp(— i%(r))u(r) proving the manifest gauge 
invariance of the boundary integral equations ([l6|), (|l8|). It can also be seen easily that 
expressions (pSj) and (29) are related by a permutation of r and Tq with subsequent 



complex conjugation (since G° is real), i.e. the operators are the adjoints of each 
other. The self-adjoint nature of (|27|) and ( |30|) follows likewise. 

The gauge independent part of the Green function, G°, has a logarithmic 
singularity at r = r . Its derivatives appearing in (E8j) - (p0[) can be expressed in 



t erm s of G° itself, at different energies v, and may be found in the appendix (A 
(A. 9). They are bounded as r — * ro. In that limit most of the quotients vanish for a 
smooth boundary, others tend to the curvature kq at the boundary point ro (defined 
to be positive for convex domains,) 



lim 



(r - r )n 



r^r (r-r ) 2 
As a consequence, all the terms in (|2 



lim 

r^r (r 



(r - r )n 



ft 

2 ' 



(31) 



(30) are integrable but for the one containing 



the (n no)/ Ar -singularity. The latter gives rise to the need for a finite part integral. 



3.2. The hyper singular integral operator 

For finite A the double-layer equation contains a hypersingular integral defined as 



N 



M = / ^-( d n /b - iA„ )(a„ /b G + iA„G)u 
f , 

- lim / - 

b 2 



= lim ^ -tf( d no/b - iA n ){d n/b G + iA„G )u. 



(32) 



We want to replace the integrand by its limiting form. To this end the boundary is 
split into the part 7 ce within an (ce)-vicinity around ro (with arbitrary constant c) 
and the remaining part r ce , 

dr 



= lim 



, 2 (d no/b - iA„ )(<9 n/6 G + iA„G )u 

Ice 



Ice 



dr 

IP 



(dno/b ~ iA„ )(9 n/6 G + iA„G )(u - u ) 



(33) 



w 



dr 

J? 



{dn /b ~ iA„ )(5 n/& G + iA„G 



with uo := u(ro). For sufficiently small e the boundary piece "f ce may be replaced by 
its tangent and the Green function by its asymptotic expression, cf. appendix A. This 
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way the third integral in 
dT 
52" 



J) may be evaluated to its contributing order, 
G + iA„cf) 



(34) 



Ice 
1 

-in 
1 

2^ 



as cos I , „ s 



V b 2 



n x r 

b 2 



0(e 2 loge) 



-s) 
1 1 



-2 



+ 4 



(s 2 + e 2 ) 2 



7r ce cr 



1 



■ 0(e loge) 



1 1 

7T C£ 



0(e 2 loge) 



0(e 2 loge). 



Here, the explicit form of the integrand was obtained from (|3C|) by the replacement 
i"o — ► r^. The last approximation in (34) holds because c may be chosen arbitrarily 
large. In a similar fashion it can be shown that the second integral in (33) is of order 
0(e). In the first integral we may replace (again for large c) the integrand by its limit 
because e is small compared to min(|r — tq|) = c e. Therefore, the limit in ( |32| ) may 
be expressed as 



AT 



-jj(d no /b - iA„ )(<9„/ 6 + iA n )Gi 



b 2 
lim 



-p-(<9n /& ~ iA„ n )(a„/ b G + iA„G)u 



f 



u - 



(35) 

LJr e "" ^e. 
where we replaced ce by e. This equation defines the finite part integral. It completes 
the derivation of the boundary integral equations and we may now turn to the question 
of how to solve them. 



4. Solving the integral equations 



As shown above, the integral equations (|16| ) and (18) may be used to compute spectra 
of magnetic billiards. However, the corresponding integral kernels are not yet suitable 
for numerical evaluation. In this section we show how their asymptotically singular 
behaviour may be separated and be treated analytically. 

In the following the combined integral equation as defined by ( |26| ) will be 
considered. The corresponding expressions for the pure double layer or single layer 
case may be obtained easily by setting a — or a -1 =0, respectively. We also take 
the opportunity to regularize the integral equations. At the energies of the Landau 
levels, f„ = n + |,ii £ No, they are defined only by the limit v — > v n , so far. This 
is because the Green function is singular at the energies v n . These simple poles are 
removed by multiplying the equations with cos(-7w) and taking the limiting values at 

Vn- 

For convenience we assume A to be constant on T and the domain T> to be simply 
connected. Let its boundary of length L — \T\ be parameterized by the arc length s, 

dr ( s ) ._ u„\ _ f- n v( s ) 



r : s G [0; L] ^ r(s) with — ^ := t(s) = "V , (36) 

as \ n x (s) J 

which allows to write the (regularized) integral kernel 

q(s,s ) := cos(7r^)[Q|j > 1 (r s ;r S0 ) +iaQ^(r s ;r 5 . ) 

TA(QS(r s ;r so )+foQN (rs;rso) )] (37) 

with r s :— r(s). After an expansion of the boundary around r(so), 

r(a) =r + (s-s )t - y (s - s ) 2 n + o((s - s ) 3 ) , (38) 
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one obtains, observing (|2 

■ r s xr 



q(s + s',s ) :-- 



b 2 



(|30"|), the asymptotic behaviour for small s' — s — sq, 
cos(niy) —1 

-2v 



p + ia T + (a - «o)^)] L -(^) 

C0S(7T^) 



k tA - 



Hp 



4tt 



+ 0(s' 2 logs' 2 ) 



(39) 



The necessary asymptotic expansions for the gauge-independent part of the Green 
function and its derivatives may be found in the appendix. L„ describes the 
asymptotically logarithmic form of the Green function and is defined in (A. 7). Note 
that due to the quotient 1/s' the expansion of zd z G° contributes up to and including 
order 0(s' logs' ). Similarly, the second order term of nno = 1 - \k 2 s' 2 + 0(s' 3 ) 
enters with the effect of cancelling another term. 

As apparent from (|3^), the singularities of the integral kernel are well described 
by the functions 

■-t-zZp-is-so) cos(tt^) -1 



m(s,s ) := TAe 1 ' 
and, for the logarithmic part, 

■ tpxro 



2tt (s-s ) 5 



(40) 



l(s,s ) 



e- b*—( s ~ S0 >L 
. (s - s ) 



)T. ( ( S ~ S rf 



b 2 



b 2 



, i ■ \ ( S Sq ) 

T A ( -p+ta-iKo)— - 2 — 



(41) 



It is important to include the terms of order 0(slog(s 2 )) to ensure that the smooth 
integral kernel defined as 

k(s, s ) := q(s, s ) - g(s - s ) [l(s, s ) + m(s, s )] (42) 

is differentiable at s = so (provided the curvature is continuous). Here, g(s') is a 
window function (with g(Q) = 1) which smoothly switches off the singular functions 
for |s'| > and vanishes beyond some small, suitably chosen a. Figure |l| shows the 
smooth as well as the original kernel for a typical choice of the boundary and the 
energy. 

The solution u(s) of the boundary integral equation is periodic and may therefore 
be expanded in a Fourier series, 



u(s) e 



E 



u e e 



(43) 



As mentioned above, we include the phase due to the gauge freedom \ which amounts 
to the choice of the symmetric gauge for the actual calculation. Within the Fourier 
representation the Fredholm determinant may be written in the form 



det 



K 



kt 



M 



kt 







(44) 



with c := (=p| — i^aA) cos(7rz/). It consists of a double Fourier integral over the smooth 
kernel, 



Km 



ds ds e 27ri (^-*ofc)/£ k(SjSo) 



t 2 



(45) 
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Figure 1. (a) Real and (c) imaginary part of the smooth combined integral kernel 
Q43) for fixed so an d the case of Neumann boundary conditions. We choose p = 0.6 
and an elliptic domain (of eccentricity 0.8 and area A = n, centered on (0.5, 0.25)) 
at v = 19 corresponding to the energy of the ~ 1000 th interior eigenstate. The 
boundary point so = is that of largest curvature. Tie magnifications (b) and 
(d) around s' = include the original singular kernel (BTj) as a dashed line. 



and two single Fourier integrals, 

L M := f ds e 2nis ^ ~ fc )/ £ L,( So ) and (46) 

M M := ^dso e 2nis °( £ ~ k )/ L M,( So ). (47) 

Here, L^(so) and M^(so) are (finite part) Fourier integrals over the asymptotic 
singularities, 

I*(*o) = f&s' e 27T[is '/ L g(s') l(s + s ) and (48) 

per 

M e (s ) =f da' e 2nils '/ L g(s') m(s + s', s ). (49) 

J —a 

They may be calculated analytically, for a suitable window <?, yielding smooth 
functions of so- In appendix B the results can be found for 

g(s') := cos 2 (|^) (e(.s' -a)- Q(s' + a)) , (50) 

where is the Heaviside step function. With this choice of the window function they 
are given in terms of elementary functions and may be evaluated easily. Having treated 
the (hyper-)singular features of the boundary integrals analytically, the remaining 
problem can be solved efficiently by numerical means. 
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Figure 2^(a) Matrix Kjy> + l^ki corresponding to the same choice of parameters 
as in Fig. hi and Dirichlet boundary conditions. Shown are the absolute values of 
the matrix elements along its diagonal and neighboring diagonals. Apart from the 
diagonal, appreciable values of the matrix + Ljy — £ c S^-i are localized within 
a sub-block which allows safe truncation. The vertical lines indicate the typical 
size after truncation, (b) The three smallest singular values of the matrix around 
u = 19 (at constant p = 0.6 corresponding to roughly the 1000 th eigenvalue.) 
The minima of the smallest singular value (solid line) determine the spectrum to 
a high accuracy. 



5. Numerical Analysis 

In the following, we describe shortly some aspects of the numerical treatment and 
discuss the question of numerical accuracy. 

The evaluation of the remaining Fourier integrals (Em) - (47) must be performed 



numerically. Since the integrands are well-behaved this may be done by dividing the 
boundary into ./V equidistant pieces and approximating the integrand at each one by 
its value at the mid-point. The summations may be performed by a Fast-Fourier 
algorithm. For large enough N this simple method is more effective than any attempt 
to evaluate the highly oscillatory integrals ( fi5| ) - ( |47| ) by more sophisticated schemes. 

Due to the Fourier representation the resulting large N x N matrix has a partly 
diagonal structure, cf. Fig. ||(a). There are off-diagonal elements of appreciable value 
only within a sub-block the size of which is independent of N. Outside of the sub-block 
essentially only the diagonal elements are occupied (the values decay rapidly as one 
leaves the diagonal.) It is the bulk wave functions which are given by the null vectors 
corresponding to the latter diagonal Fourier components. These components do not 
contribute to the other states since they are not coupled to them. As a consequence, 
the restriction of the matrix to the above-mentioned sub-block at most removes bulk 
states, if they exist, out of the numerically calculated spectrum without affecting other 
states. Generically, one is not particularly interested in these states whose energies are 
exponentially close to the Landau levels. Since the spectrum is modified at most in a 
well-controlled way, it is permissible to truncate the matrix to a smaller size Nt Iunc . 

A small complication arises in the case of finite A. Due to the hypersingular 
part the diagonal Fourier elements increase linearly as \£\ — * oo, cf. (|B.5|). The 
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above statements apply in this case after dividing the matrix ( |44| ) column-wise by 
the asymptotic factor 



•<t x r > 1,2 , /Si(7r) 

Z7T- 



1/2 



62 .-^ ■ v a ,j. (51) 

Here, < to xro > is the average (the th Fourier component) of the function t(so) xr(so) 
defined on the boundary. 

The calculation of the spectrum amounts to finding (all) the zeros of the complex- 
valued determinant ([ii] ) in a given energy range. Numerically, this is the most 
expensive task, scaling as iV trunc . Since the computation of the determinant tends 
to be unstable around its zeros it is more advantageous to employ a singular-value 
decomposition of the matrix which is stable in any case. The vanishing of a singular 
value indicates a defective rank of its matrix. Due to roundoff errors these non- negative 
quantities are always greater than zero. However, the spectral points are very well 
defined by the sharp minima of the lowest singular value as a function of v, cf. Figure 
||(b). The detection of near degeneracies is made appreciably easier if one monitors 
the next smallest singular values, too. 

In order to calculate the wave function tf>Q = ip(rQ (fc T) away from the boundary 
one may use directly equation (|l2|). The gauge invariant gradient of the wave function, 
7 := Vr/bV'o — iAo"0o needed for the current density J = Im^gTo] * s obtained from 
the same equation after the application of the operator v ro /j — iAo. 

^o = ±y^[±£(d n/b G + iA„G)-G]u, (52) 
— [ ± £ (v ro/6 - iA )(a n/6 G + iA„G) - (v ro/b G - iA G)] w (53) 

Since the integrands are not singular for To ^ T the integrals may be approximated 
by a discrete sum over the N boundary elements without further ado. The densities 
of other observables can be obtained by similar boundary integrals. 



5.1. Convergence and Accuracy 

Careful numerical tests indicate that the precision of the calculated spectra and wave 
functions is determined almost exclusively by the dimension N of the initial matrix. 
In Figure |^(a) we show how the energies converge exponentially as N increases. At the 
same time, the calculated spectra are found to be numerically invariant with respect 
to other parameters such as a, a, iVt runc , and in particular the location of the origin. 

Reasonable choices for a and a are a = vj (2p) and a = b. The location and size of 
the sub-block is best determined in terms of an averaged column norm. The resulting 
spectra are independent of -ZV trunc provided it exceeds a critical value. Here, the 
position of the origin is relevant, because the calculation of the spectral determinant 
(Q), in particular its analytical parts, must be performed in a specific gauge. The 
choice in favor of the symmetric gauge is made in (Q) where the remaining gauge 
freedom \ is absorbed into the solution vector. As a consequence of the resulting 
distinction of the origin, the spectral determinant is no longer translationally invariant. 

As a result, the size of the truncated matrix depends on the choice of the origin. 
For example, the values in Fig. ||(a) belong to an ellipse centered at the origin. With 
an ellipse displaced by (2,1) one obtains the same relative errors for N = 600 . . . 2400 
(not shown, one would not see a difference) with truncation sizes larger by 50%. In 
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Figure 3. Errors of the ~ 1000 th interior eigenvalue at p = 0.6 as a function 
of the boundary discretization N. (a) Approximate relative error for the elliptic 
domain of Fig. ^](b) (the Dirichlet state closest to v = 19) Here, the energy for 
N = 2600 was taken as reference. The numbers indicate the matrix dimension 
after truncation which determines the numerical effort. They increase only weakly 
with N. (b) Exact relative error of the exterior Neumann energies of a typical 
edge state ( A, V) and a typical bulk state (□, ()) as a function of N. Here, 
we use a circular domain (of area A = it) which allows to determine the exact 
energies (f e dge — 19.0294509, i^buik — 19.4816851) independently. The center of 
the domain is placed at the origin (A, □) and at the point (3,0) (V, 0)i respectively. 
One observes that the displacement does not affect the error of the edge state but 
increases the error of the bulk state energy systematically. [Note, that the graphs 
do not have the same scale. 1 



order to minimize the numerical effort it is therefore advantageous to choose the origin 
in the center of the domain considered. 

The fact that bulk states are more strongly affected by the truncation is seen 
in Fig. ||(b) where exterior Neumann states of a circular domain are compared after 
displacement by 3 radii. Since the disk is a separable problem, we can check here 
against the exact energies (obtained as the roots of an analytical function.) Note, 
that the calculation of the hypersingular integral introduces no additional error. 

The only precise published calculations for a nontrivial shape known to us are 
the results of Tiago et al who give the first twenty Dirichlet levels for an ellipse of 
eccentricity 0.8 and area A = tt at constant b 2 = 2/25 (missing one symmetry class.) 
Our method is able to confirm their results to all given seven digits (apart from 
occasional differences in the last digit by one). For reference, we note the energy of 
the approximately thousandth state (the one closest to v = 80) which we calculate to 
be v ~ 79.9362(6). The expected error is less than 0.1% of the mean level spacing. 

6. Numerical Results 

In the following we demonstrate the performance of the described method by exhibiting 
some numerical results on magnetic billiards which have been inaccessible by other 
methods. 
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Figure 5. Fluctuating part of the spectral staircase, Nfl uc := N(f) — N for the 
asymmetric stadium at p = 1.2. The displayed range contains the first 5000 points 
in the spectrum. The heavy line is a running average over 250 neighboring points. 
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6.1. Spectral statistics 

We start by considering spectral statistics based on large data sets of calculated 
spectral points. As explained in the introduction, the spectra are defined in the 
semiclassical direction 6 — >■ keeping the cyclotron radius p constant. This way the 
underlying classical dynamics is fixed. For classically hyperbolic systems one expects 
Random Matrix Theory (RMT) to reproduce the spectral statistics. 

We use the two domains described in Fig. [|. One is an asymmetric version of the 
Bunimochwich stadium billiard (r\ = 0.75, r2 = 0.25, .A = 5.39724). In the magnetic 
field its dynamics is free of unitary symmetries but contains an anti-unitary one (time 
reversal and reflection at y = 0.) On the other hand, the skittle shape (made up of the 
arcs of four symmetrically touching circles, r x = 1.0, = 0.5, A = 4.33969) does not 
have any symmetry. We choose it because it generates hyperbolic classical motion even 
for small cyclotron radii p > 1, according to a recent theorem (6). The asymmetric 
stadium could not be proven to be hyperbolic although we find no numerical evidence 
for systematic deviations from the RMT behavior (see below.) 

We calculated 5300 and 7300 consecutive interior Dirichlet eigenvalues at p = 1.2 
for the asymmetric stadium and the skittle shaped domain, respectively. It should be 
noted, that states of much higher ordinal number can be computed at little cost with 
the present method. The time consuming task is rather to find all energies Vi = p 2 /bf, 
including the near-degenerate ones, in a given interval. 

A quantity which sensitively indicates whether spectral points were missed is the 
fluctuating part, Nfl uct , of the spectral staircase function 

N(i/) := V 6(i/ - Vi ) = N(i/) + N fluct (^). (54) 
* — *i 

As shown recently ^7j, its smooth part coincides with the non- magnetic one in its 
leading terms. In our units and for Dirichlet boundary conditions they read 

ffM^^-r"^' (55) 

p z 7r Zirp o 

where A is the domain area and L the boundary length. The constant, which contains 
an integral over the boundary curvature, is the same for the shapes considered. Figure 
H displays the fluctuating part of the number function Nfl uc := N(f) — N for the 
asymmetric stadium. It proves that the spectrum is complete. A similar result is 
obtained for the skittle shaped domain (not shown.) 

The large spectral intervals at hand allow us to calculate directly some of the 
popular spectral functions. Due to the underlying classical chaos and the symmetry 
properties mentioned above one expects the statistics of the Gaufiian Orthogonal 
Ensemble (GOE) for the asymmetric stadium and the Gaufiian Unitary Ensemble 
(GUE) for the skittle. Figure || shows the distributions of nearest neighbors P(s) 
of the unfolded spectra. Indeed, one finds excellent agreement with Random Matrix 
Theory. The differences between the numerical and the RMT cumulative functions 
J(a) = f° P(s')ds' stay below 2%. 

A function which characterizes the spectrum much more sensitively than P(s) is 
the form factor K(t), i.e. the (spectrally averaged) Fourier transform of the two-point 
autocorrelation function of the spectral density |2?| . Figure ^ gives the spectral 
form factor together with the RMT results. We find very good agreement. One would 
expect systematic deviations at small r due to the contributions of single short periodic 
orbits. These cannot be resolved with the present size of the spectral interval, though. 
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Figure 6. Nearest neighbor distributions of the asymmetric stadium (left) and 
the skittle shaped domain (right), at p = 1.2. The histograms should be compared 
to GOE and GUE predictions of Random Matrix Theory, respectively (heavy 
lines.) The monotonic lines give the corresponding cumulative quantities. Their 
differences are reported in the insets. 



Since most other popular spectral measures like Dyson's A3 statistic are functions of 
the form factor there is no need to present them here. 

The good agreement with RMT is not only a consequence of the large spectral 
intervals the statistics are based on. It is equally important that the spectra are 
defined at fixed classical dynamics. Had we calculated the spectra at fixed field, 
they would have been based on a classical phase space that transforms from a near- 
integrable, time-invariance-broken structure to a hyperbolic time-invariant one as p 
increases with energy. This transformation of spectral statistics from GOE to GUE 
as the field is increased was studied in [[l3| [l4|, |l5) . 

6.2. Wave junctions 

We proceed to present a selection of wave functions calculated in the semiclassical 
regime. We start with those of the skittle shaped domain choosing again p = 1.2. 
This ensures that the corresponding classical skipping motion is hyperbolic in the 
interior, as well as in the exterior. 

Figure ||(a) shows the density plot of a typical interior wave function around the 
one-thousandth eigenstate. As expected for a classically ergodic system, it spreads 
out throughout the whole domain but is not completely featureless. Occasionally, 
one may also find bouncing-ball modes, i.e. wave functions localized on a manifold of 
marginally stable periodic orbits. One such wave function is given in Fig. ^|(b). It 
belongs to a family of 2-orbits. 

A typical exterior wave function with an energy close to that of Fig. ^|(a) is 
displayed in the middle row of Figure ||, at the same scale (c) and a larger scale (d). 
One observes that in the vicinity of the boundary it behaves quite similar to an interior 
function. On a larger scale, the wave function decays after a distance smaller than two 
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Figure 7. Spectral form factor of the asymmetric stadium (left) and the skittle 
shaped domain (right), based on 5300 and 7300 spectral points, respectively. The 
heavy lines display the same data after stronger spectral averaging. The random 
matrix result for the GauBian Orthogonal and the Gaufiian Unitary Ensemble, 
respectively, is indicated by the dashed lines. The insets show the regions of small 



cyclotron radii. In this region circular structures are faintly visible with the radius of 
the classical cyclotron motion. 

The bottom row of Figure |^ shows a quite different exterior state with an energy 
close to that of a Landau level. It is a bulk state. A typical feature is the fact there 
are no large amplitudes close to the boundary. Rather, one finds a ring of increased 
amplitude encircling the domain. Another ring surrounds the domain at a distance 
of 2p. This double-ring structure moves outwards as one goes through the series of 
states with energies increasingly close to the Landau levels. Semiclassically, it can be 
understood as being made up of a superposition of cyclotron orbits. This becomes 
even more clear in the following where we consider a more symmetric shape of the 
boundary. 

For the second set of wave functions we choose an elliptic domain (of eccentricity 
0.8 and area 7r) at a small cyclotron radius p = 0.6. The classical dynamics is mixed 
chaotic in this case || . Going to the extreme semiclassical limit - the ten-thousandth 
interior eigenstate - we expect the wave functions to mimic the structures of the 
underlying classical phase space. 

Indeed, Figure ^(a) displays a wave function which is localized along a stable 
interior 6 x 6-orbit. Note that the wave nature of the eigenstate is still visible at points 
where the trajectory crosses with itself, in particular at the shallow intersections close 
to the center. 

Since p is small enough to allow closed cyclotron orbits within the ellipse, we find 
bulk states also in the interior, see Fig. ^(b) for an example. Again it is semiclassically 
described by a superposition of closed cyclotron orbits. This can be seen clearly from 
the current distributions which are given in the bottom row of Fig. for the edge state 
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Figure 8. Interior and exterior wave functions of the skittle shape at p = 1.2 
around the one-thousandth interior eigenstate. The plotted shade is proportional 
to the thin lines indicate the boundary T. Energies: (a) v ~ 32.98804, (b) 
v ~ 33.12033, (c,d) v ~ 32.84740, (e,f) y ~ 32.50073. 




Figure 9. Wave functions (a,b) and current distributions (c,d) in an elliptic 
domain at p = 0.6, around the ten-thousandth interior eigenstate, with energies 
v ~ 60.06026 (a,c) and v ~ 60.50030 (b,d). 



(c) and the bulk state (d), respectively. Here, the length of the arrows is proportional 
to the amplitude of the current density. 

Similar semiclassical states can also be found in the exterior, as displayed in 
Figure [Io|. The edge state, Fig. |T^(a), obviously belongs to an 8 x 6-orbit. Like all 
edge states it is distinguished from a typical bulk state, cf. Fig. |ic](b), by the finite 
current it carries around the domain. In contrast, the bulk state with its counter- 
running current densities has no net current along the boundary, cf. Fig. |Io|(c) and 
Fig. 0(d). 

We emphasize that all the wave functions and current distributions shown above 




Figure 10. Exterior wave functions (aj)) and current distributions (c,d) at 
p = 0.6 and at similar energies as in Fig. H, v ~ 60.13634 (a,c) and v ~ 60.50049 
(b,d). 
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Figure 11. Parametric dependence of the exterior spectrum on the boundary 
condition (for the asymmetric stadium at fixed b = 0.25). The parameter 
A interpolates between Dirichlet (A = 0) and Neumann (A = 1) boundary 
conditions. The right graph shows details around the forth Landau level. 

are calculated throughout the entire displayed area. They turn out to be numerically 
zero in the complementary domains as expected from the theory Consequently, the 
type of a solution of a single integral equation can be inferred by calculating the wave 
function. 

6.3. General boundary conditions 

So far, only Dirichlet boundary conditions were considered. As a last point we show 
the parametric dependence of a spectrum on the type of boundary conditions. Figure 
|lT| presents the exterior spectrum of the asymmetric stadium, calculated at fixed 
b = 0.25. The value A is taken constant along the boundary and parameterized by a 
number Ae [0; 1], 

A = -^-tan(fA). (56) 

Here A is chosen negative to ensure that the transformation from Dirichlet (A = 0) to 
Neumann (A = 1) boundary conditions is continuous. For positive A this would not 
be the case, which is a restriction similar to the one for the field free case (2^] . 

The energies clustering around the Landau levels v n = n + i , n e No belong 
to bulk states. One observes that they are lifted from the Landau levels into higher 
energies at Dirichlet boundary conditions, whereas in the Neumann case they are 
shifted to smaller energies. A semiclassical theory which describes the exponential 
approach of the bulk states to the Landau levels and its transition as a function of A 
will be published elsewhere. 
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Conclusions 

The main theoretical result presented here is the finding that the two boundary integral 
equations of the billiard problem admit spurious solutions in the magnetic case, and 
how those are identified and removed. 

An important implication concerns the semiclassical theory of magnetic billiards. 
The trace formula for the semiclassical quantization of the field free case is based on a 
(double layer) boundary integral equation [ |30| . If one tries to repeat the derivation for 
the magnetic case, problems should arise since the equation does not give a sufficient 
condition. Indeed, a (Balian-Bloch-like) derivation of the trace formula for magnetic 
billiards does not exist. We emphasize that the starting point should be the gauge- 
invariant formulation of the boundary integral equation, as presented above. Only 
then are the spurious solutions gauge invariant, have a physical interpretation, and 
can be taken into account systematically. 

We have shown how a precise and efficient computational method for the 
calculation of spectra and wave functions can be based on a combination of boundary 
integral equations. It allows to obtain the exact spectra and wave functions at energies 
and fields inaccessible hitherto. 

The possibility of calculating the exterior spectra as well, brings up the question 
on how interior and exterior spectra are related. Here, a problem is the existence 
of an infinity of bulk states which do not have much physical relevance but prevent 
the spectral number function from being well-defined. Our calculation of the exterior 
level dynamics as a function of the boundary condition shows that the edge states 
and bulk states differ in their sensitivity to the boundary condition. In a forthcoming 
communication we will propose a definition of the spectral edge state density based 
on this observation. 
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Appendix A. The free magnetic Green function 

The free magnetic Green function was derived in |3l[ fL7| by angular momentum 
decomposition. Here, we show how it is obtained by directly performing the Fourier 
transform of the time evolution operator |32| 

U(r, r ;*) = -±4~^ e J (^ + - xM) (A . 1} 

27ri b z sm(wi) 

which yields both, the gauge dependent and the gauge independent part in a 
straightforward manner. We have to evaluate 



G(r;r ) = h l J\ut) e iEt / h V(r,r ;t) 



- X(r) + xM) G ((£^0)!) (A . 2) 
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assuming that the energy v has a small imaginary part to ensure convergence. For 
the gauge independent part one obtains 
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(A.3) 



Air ' y2 

where we used a logarithmic representation of the inverse cotangens and the reflection 



relation r(| — v)T{^ + v) cos(-7w) = tt. The last equality in (A.3) holds since the 
expression in square brackets may be deformed to the (complex conjugate of the) 
contour integral found in |3^, eq. (5.1.6)]. It gives the (real valued) irregular Whittaker 



function W 34 1 (multiplied by z 2 ) in an integral representation that is valid even 
for positive v. The regularized version of G , 



G°(z) := lim cos(7r/z)G°(z), 



f_L — >V 



(A.4) 



reads in terms of the more common irregular confluent hypergeometric function U 



-z/2 



4tt r(u+i) 



m-v,l;z). 



At small distances, it has the logarithmic form, 

G°(z) = K(z) + O(zlogz) asz^O, 
cos(7rf) 



K(z) 



Air 



-(log(z) + *(i + i/)-2*(l))- 



sin(7rt/) 



(A.5) 

(A.6) 
(A.7) 



where VP is the Digamma function. 



A.l The derivatives and their asymptotic behavior 

Employing the differential properties of the confluent hypergeometric function |34j we 
can express the derivatives of by the function G^ itself. 

z-^GUz) = - (i - «/)(G° + G°_ x ) - |G° (A.8) 
^ 2 ^G«(z) = (| - - ^)(G° + 2G°_ X + G°_ 2 ) 



+ z(I-z,)(G°+G°_ 1 ) + ^ r G° 



(A.9) 



In section [l| we need their asymptotic expansions, 

z AqO(z) = C ° S f 7r ' y) [1 - z v ( log(z) + *(i - v) - 2* (1) - 1)1 + 0(z 2 log z) (A.10) 
dz 47T L v z /j 



dz ; 



:G°(z) = - 



Air 



O(zlogz) as z 



0. 



(A.ll) 
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These were deduced from the logarithmic representation of U in terms of the regular 
Kummer function pj, eq. (13.6.1)]. 



Appendix B. Analytical calculation of the singular integrals 



The Fourier integrals fl48|), (49) depend on the the window function g. Our choice 
is (^0|) which switches off the asymptotically singular functions m and 1 sufficiently 
smoothly. For the logarithmic integrals one obtains 
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where f2;(so) = 2n£/L + to fc * r ° , (p = Qi(sq) a, ip^ 1 — tp± tt, and Si is the Sine Integral. 
The finite part integral reads 
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Asymptotically, 

M,( S o)^TA^^r»/( S o)sgn(£) as \£\ -» oo. 
Note that with choice ([50]) the limit of the remaining kernel is 
q(s, so) - - s ) (l(s, s ) + m(s, s )) 
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which is noi just the constant part of (p9|) but contains a term which depends on a . 
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Figure 12. (a) Gauge independent part of the regularized Green function at 
v = 57.75. It has a logarithmic singularity at r = and decays exponentially 
for r > 2p. (b) In the transition regions between oscillatory, transient, and 
decaying regimes the asymptotic expressions to third order are not valid (chain, 
dotted, dashed line respectively.) (c) Here, on may interpolate using an uniform 
approximations to the irregular Whittaker function (solid line.) 



Appendix C. Numerical evaluation of the Green function 

We are not aware of any published numerical procedure to evaluate the irregular 
confluent hypergeometric function U if both, the (energy) parameter and the variable 
are large. It seems that presently only the Mathematica software (Wolfram Research 
Inc.) is able to compute the function, at least for moderately large v. Even this 
sophisticated system fails for v > 75. Anyhow, it is not an option to use it for serious 
numerical calculations since the evaluation takes a prohibitively long time. 

Therefore, wc describe our method to compute the gauge independent part of 
the regular Green function in more detail. For low energies v < 12, the function 
U(l/2 — v, 1; z) may be easily calculated by its series representation J34|, eq. (13.1.6)]. 
i.e. in terms of the regular confluent hypergeometric function iFi. For very large z an 
asymptotic expansion in terms of 2F3 may be employed |35[ eq. (6.7.1)]. 

For energies v > 12 the numerical convergence of the series expression deteriorates 
strongly in some intervals of the z range (starting at z w 2i>) . Fortunately, a number 
of rather complicated asymptotic expansions for the irregular Whittaker function 
exist eqs. (8.1.5), (8.1.10), (8.1.18a)] which are to third order in the large 
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parameter v. These expressions are based on saddle point approximations of a defining 
integral. Together with eq. (13.5.15)] they correspond to the changing logarithmic, 
oscillatory, transient, and exponentially decaying behavior of the Green function as the 
distance z increases. For most values of z they allow to calculate the Green function to 
a reasonably high precision and with acceptable numerical effort. However, between 
the ranges of validity of the different asymptotic expressions there are small gaps where 
no formula is appropriate, cf. Figure [l^. In the gap between the logarithmic and the 
oscillatory domains, which is at small z, one may employ the series summation even 
for large v ^> 12. For the two gaps between the oscillatory, the transient, and the 
exponential regimes, which are around z ~ Av this is possible only up to, say v = 16. 
For larger v we interpolate between adjacent regions of validity employing the uniform 
approximation of the irregular Whittaker function around the classical turning point. 
Neglecting higher orders in v, the resulting expression for the Green function reads, 

G v (z)*C — { \ q)6 Ai(sgn(z-z )(|g) § ) (C.l) 
\z z — Avz — 1| * v ' 

where Ai is the regular Airy function and 

w JJ 2 5 Vz l + 2vz 



(IT fz—Zv\\ 1, /Zq 1 + Ivz + w\ 1 

v\ atan M — log w it z < zn 

V2 V to J J 2 B V z l + 2z/z / 2 ° 

1 1 /2vZ + l\ 7T (z — 2v + w\ 

-w + -atan( I ---i/ log I if z > z Q 

2 2 V w / 4 \ zq — 2v / 



(C.2) 



with z = Av (i + Wl + jir) and to = y/\z 2 -Avz- 1|. (C.3) 

The constant C may be calculated for values of z where the saddle point expressions 
are valid and is interpolated linearly within the gaps. 

The thresholds mentioned above are a reasonable compromise between cost and 
precision. We observe a peak numerical error (minimum of relative and absolute) of 
6.5 x 10 -5 at v — 22 by comparison with the results of Mathematica which are assumed 
to be exact for v < 70. For increasing v the numerical error decreases monotonically 
which allows us to estimate it to smaller than 3.7 x 1CP 5 for v > 70. It was checked 
that numerical errors of that order do not affect the results shown in section 6. 
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